home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / tos / g__~1 / gplibs15.zoo / iodtoa.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1992-07-23  |  74.0 KB  |  2,405 lines
  1. #include <ioprivat.h>
  2. #ifdef USE_DTOA
  3. /****************************************************************
  4.  *
  5.  * The author of this software is David M. Gay.
  6.  *
  7.  * Copyright (c) 1991 by AT&T.
  8.  *
  9.  * Permission to use, copy, modify, and distribute this software for any
  10.  * purpose without fee is hereby granted, provided that this entire notice
  11.  * is included in all copies of any software which is or includes a copy
  12.  * or modification of this software and in all copies of the supporting
  13.  * documentation for such software.
  14.  *
  15.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  16.  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  17.  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  18.  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  19.  *
  20.  ***************************************************************/
  21.  
  22. /* Please send bug reports to
  23.         David M. Gay
  24.         AT&T Bell Laboratories, Room 2C-463
  25.         600 Mountain Avenue
  26.         Murray Hill, NJ 07974-2070
  27.         U.S.A.
  28.         dmg@research.att.com or research!dmg
  29.  */
  30.  
  31. /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  32.  *
  33.  * This strtod returns a nearest machine number to the input decimal
  34.  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  35.  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  36.  * biased rounding (add half and chop).
  37.  *
  38.  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  39.  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  40.  *
  41.  * Modifications:
  42.  *
  43.  *      1. We only require IEEE, IBM, or VAX double-precision
  44.  *              arithmetic (not IEEE double-extended).
  45.  *      2. We get by with floating-point arithmetic in a case that
  46.  *              Clinger missed -- when we're computing d * 10^n
  47.  *              for a small integer d and the integer n is not too
  48.  *              much larger than 22 (the maximum integer k for which
  49.  *              we can represent 10^k exactly), we may be able to
  50.  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  51.  *      3. Rather than a bit-at-a-time adjustment of the binary
  52.  *              result in the hard case, we use floating-point
  53.  *              arithmetic to determine the adjustment to within
  54.  *              one bit; only in really hard cases do we need to
  55.  *              compute a second residual.
  56.  *      4. Because of 3., we don't need a large table of powers of 10
  57.  *              for ten-to-e (just some small tables, e.g. of 10^k
  58.  *              for 0 <= k <= 22).
  59.  */
  60.  
  61. /*
  62.  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  63.  *      significant byte has the lowest address.
  64.  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  65.  *      significant byte has the lowest address.
  66.  * #define Sudden_Underflow for IEEE-format machines without gradual
  67.  *      underflow (i.e., that flush to zero on underflow).
  68.  * #define IBM for IBM mainframe-style floating-point arithmetic.
  69.  * #define VAX for VAX-style floating-point arithmetic.
  70.  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
  71.  * #define No_leftright to omit left-right logic in fast floating-point
  72.  *      computation of dtoa.
  73.  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
  74.  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  75.  *      that use extended-precision instructions to compute rounded
  76.  *      products and quotients) with IBM.
  77.  * #define ROUND_BIASED for IEEE-format with biased rounding.
  78.  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  79.  *      products but inaccurate quotients, e.g., for Intel i860.
  80.  * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
  81.  *      integer arithmetic.  Whether this speeds things up or slows things
  82.  *      down depends on the machine and the number being converted.
  83.  * #define KR_headers for old-style C function headers.
  84.  */
  85.  
  86. #ifdef DEBUG
  87. #include <stdio.h>
  88. #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  89. #endif
  90.  
  91. #include <stdlib.h>
  92. #include <string.h>
  93. #define CONST const
  94.  
  95. #include <errno.h>
  96. #include <float.h>
  97. #ifndef __MATH_H__
  98. #include <math.h>
  99. #endif
  100.  
  101. #ifdef atarist
  102. #define IEEE_MC68k
  103. #endif
  104.  
  105. #ifdef Unsigned_Shifts
  106. #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
  107. #else
  108. #define Sign_Extend(a,b) /*no-op*/
  109. #endif
  110.  
  111. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  112.  
  113. #if FLT_RADIX==16
  114. #define IBM
  115. #elif DBL_MANT_DIG==56
  116. #define VAX
  117. #elif DBL_MANT_DIG==53 && DBL_MAX_10_EXP==308
  118. #define IEEE_Unknown
  119. #else
  120. Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  121. #endif
  122. #endif
  123.  
  124. #ifdef IEEE_8087
  125. #define HIWORD 1
  126. #define LOWORD 0
  127. #define TEST_ENDIANNESS  /* nothing */
  128. #elif defined(IEEE_MC68k)
  129. #define HIWORD 0
  130. #define LOWORD 1
  131. #define TEST_ENDIANNESS  /* nothing */
  132. #else
  133. static int HIWORD = -1, LOWORD;
  134. static void test_endianness()
  135. {
  136.     union doubleword {
  137.     double d;
  138.     unsigned long u[2];
  139.     } dw;
  140.     dw.d = 10;
  141.     if (dw.u[0] != 0) /* big-endian */
  142.     HIWORD=0, LOWORD=1;
  143.     else
  144.     HIWORD=1, LOWORD=0;
  145. }
  146. #define TEST_ENDIANNESS  if (HIWORD<0) test_endianness();
  147. #endif
  148. #define word0(x) ((unsigned long *)&x)[HIWORD]
  149. #define word1(x) ((unsigned long *)&x)[LOWORD]
  150.  
  151. /* The following definition of Storeinc is appropriate for MIPS processors. */
  152. #if defined(IEEE_8087) + defined(VAX)
  153. #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
  154. ((unsigned short *)a)[0] = (unsigned short)c, a++)
  155. #elif defined(IEEE_MC68k)
  156. #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
  157. ((unsigned short *)a)[1] = (unsigned short)c, a++)
  158. #else
  159. #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
  160. #endif
  161.  
  162. /* #define P DBL_MANT_DIG */
  163. /* Ten_pmax = floor(P*log(2)/log(5)) */
  164. /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  165. /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  166. /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  167.  
  168. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_Unknown)
  169. #define Exp_shift  20
  170. #define Exp_shift1 20
  171. #define Exp_msk1    0x100000
  172. #define Exp_msk11   0x100000
  173. #define Exp_mask  0x7ff00000
  174. #define P 53
  175. #define Bias 1023
  176. #define IEEE_Arith
  177. #define Emin (-1022)
  178. #define Exp_1  0x3ff00000
  179. #define Exp_11 0x3ff00000
  180. #define Ebits 11
  181. #define Frac_mask  0xfffff
  182. #define Frac_mask1 0xfffff
  183. #define Ten_pmax 22
  184. #define Bletch 0x10
  185. #define Bndry_mask  0xfffff
  186. #define Bndry_mask1 0xfffff
  187. #define LSB 1
  188. #define Sign_bit 0x80000000
  189. #define Log2P 1
  190. #define Tiny0 0
  191. #define Tiny1 1
  192. #define Quick_max 14
  193. #define Int_max 14
  194. #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
  195. #else
  196. #undef  Sudden_Underflow
  197. #define Sudden_Underflow
  198. #ifdef IBM
  199. #define Exp_shift  24
  200. #define Exp_shift1 24
  201. #define Exp_msk1   0x1000000
  202. #define Exp_msk11  0x1000000
  203. #define Exp_mask  0x7f000000
  204. #define P 14
  205. #define Bias 65
  206. #define Exp_1  0x41000000
  207. #define Exp_11 0x41000000
  208. #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
  209. #define Frac_mask  0xffffff
  210. #define Frac_mask1 0xffffff
  211. #define Bletch 4
  212. #define Ten_pmax 22
  213. #define Bndry_mask  0xefffff
  214. #define Bndry_mask1 0xffffff
  215. #define LSB 1
  216. #define Sign_bit 0x80000000
  217. #define Log2P 4
  218. #define Tiny0 0x100000
  219. #define Tiny1 0
  220. #define Quick_max 14
  221. #define Int_max 15
  222. #else /* VAX */
  223. #define Exp_shift  23
  224. #define Exp_shift1 7
  225. #define Exp_msk1    0x80
  226. #define Exp_msk11   0x800000
  227. #define Exp_mask  0x7f80
  228. #define P 56
  229. #define Bias 129
  230. #define Exp_1  0x40800000
  231. #define Exp_11 0x4080
  232. #define Ebits 8
  233. #define Frac_mask  0x7fffff
  234. #define Frac_mask1 0xffff007f
  235. #define Ten_pmax 24
  236. #define Bletch 2
  237. #define Bndry_mask  0xffff007f
  238. #define Bndry_mask1 0xffff007f
  239. #define LSB 0x10000
  240. #define Sign_bit 0x8000
  241. #define Log2P 1
  242. #define Tiny0 0x80
  243. #define Tiny1 0
  244. #define Quick_max 15
  245. #define Int_max 15
  246. #endif
  247. #endif
  248.  
  249. #ifndef IEEE_Arith
  250. #define ROUND_BIASED
  251. #endif
  252.  
  253. #ifdef RND_PRODQUOT
  254. #define rounded_product(a,b) a = rnd_prod(a, b)
  255. #define rounded_quotient(a,b) a = rnd_quot(a, b)
  256. extern double rnd_prod(double, double), rnd_quot(double, double);
  257. #else
  258. #define rounded_product(a,b) a *= b
  259. #define rounded_quotient(a,b) a /= b
  260. #endif
  261.  
  262. #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  263. #define Big1 0xffffffff
  264.  
  265. #ifndef Just_16
  266. /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
  267.  * This makes some inner loops simpler and sometimes saves work
  268.  * during multiplications, but it often seems to make things slightly
  269.  * slower.  Hence the default is now to store 32 bits per long.
  270.  */
  271. #ifndef Pack_32
  272. #define Pack_32
  273. #endif
  274. #endif
  275.  
  276. #define Kmax 15
  277.  
  278. extern "C" double strtod(const char *s00, char **se);
  279. extern "C" char *dtoa(double d, int mode, int ndigits,
  280.                         int *decpt, int *sign, char **rve);
  281.  
  282.  struct
  283. Bigint {
  284.         struct Bigint *next;
  285.         int k, maxwds, sign, wds;
  286.         unsigned long x[1];
  287.         };
  288.  
  289.  typedef struct Bigint Bigint;
  290.  
  291.  static Bigint *freelist[Kmax+1];
  292.  
  293.  static Bigint *
  294. Balloc
  295. #ifdef KR_headers
  296.         (k) int k;
  297. #else
  298.         (int k)
  299. #endif
  300. {
  301.         int x;
  302.         Bigint *rv;
  303.  
  304.         if (rv = freelist[k]) {
  305.                 freelist[k] = rv->next;
  306.                 }
  307.         else {
  308.                 x = 1 << k;
  309.                 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
  310.                 rv->k = k;
  311.                 rv->maxwds = x;
  312.                 }
  313.         rv->sign = rv->wds = 0;
  314.         return rv;
  315.         }
  316.  
  317.  static void
  318. Bfree
  319. #ifdef KR_headers
  320.         (v) Bigint *v;
  321. #else
  322.         (Bigint *v)
  323. #endif
  324. {
  325.         if (v) {
  326.                 v->next = freelist[v->k];
  327.                 freelist[v->k] = v;
  328.                 }
  329.         }
  330.  
  331. #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
  332. y->wds*sizeof(long) + 2*sizeof(int))
  333.  
  334.  static Bigint *
  335. multadd
  336. #ifdef KR_headers
  337.         (b, m, a) Bigint *b; int m, a;
  338. #else
  339.         (Bigint *b, int m, int a)       /* multiply by m and add a */
  340. #endif
  341. {
  342.         int i, wds;
  343.         unsigned long *x, y;
  344. #ifdef Pack_32
  345.         unsigned long xi, z;
  346. #endif
  347.         Bigint *b1;
  348.  
  349.         wds = b->wds;
  350.         x = b->x;
  351.         i = 0;
  352.         do {
  353. #ifdef Pack_32
  354.                 xi = *x;
  355.                 y = (xi & 0xffff) * m + a;
  356.                 z = (xi >> 16) * m + (y >> 16);
  357.                 a = (int)(z >> 16);
  358.                 *x++ = (z << 16) + (y & 0xffff);
  359. #else
  360.                 y = *x * m + a;
  361.                 a = (int)(y >> 16);
  362.                 *x++ = y & 0xffff;
  363. #endif
  364.                 }
  365.                 while(++i < wds);
  366.         if (a) {
  367.                 if (wds >= b->maxwds) {
  368.                         b1 = Balloc(b->k+1);
  369.                         Bcopy(b1, b);
  370.                         Bfree(b);
  371.                         b = b1;
  372.                         }
  373.                 b->x[wds++] = a;
  374.                 b->wds = wds;
  375.                 }
  376.         return b;
  377.         }
  378.  
  379.  static Bigint *
  380. s2b
  381. #ifdef KR_headers
  382.         (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
  383. #else
  384.         (CONST char *s, int nd0, int nd, unsigned long y9)
  385. #endif
  386. {
  387.         Bigint *b;
  388.         int i, k;
  389.         long x, y;
  390.  
  391.         x = (nd + 8) / 9;
  392.         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
  393. #ifdef Pack_32
  394.         b = Balloc(k);
  395.         b->x[0] = y9;
  396.         b->wds = 1;
  397. #else
  398.         b = Balloc(k+1);
  399.         b->x[0] = y9 & 0xffff;
  400.         b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
  401. #endif
  402.  
  403.         i = 9;
  404.         if (9 < nd0) {
  405.                 s += 9;
  406.                 do b = multadd(b, 10, *s++ - '0');
  407.                         while(++i < nd0);
  408.                 s++;
  409.                 }
  410.         else
  411.                 s += 10;
  412.         for(; i < nd; i++)
  413.                 b = multadd(b, 10, *s++ - '0');
  414.         return b;
  415.         }
  416.  
  417.  static int
  418. hi0bits
  419. #ifdef KR_headers
  420.         (x) register unsigned long x;
  421. #else
  422.         (register unsigned long x)
  423. #endif
  424. {
  425.         register int k = 0;
  426.  
  427.         if (!(x & 0xffff0000)) {
  428.                 k = 16;
  429.                 x <<= 16;
  430.                 }
  431.         if (!(x & 0xff000000)) {
  432.                 k += 8;
  433.                 x <<= 8;
  434.                 }
  435.         if (!(x & 0xf0000000)) {
  436.                 k += 4;
  437.                 x <<= 4;
  438.                 }
  439.         if (!(x & 0xc0000000)) {
  440.                 k += 2;
  441.                 x <<= 2;
  442.                 }
  443.         if (!(x & 0x80000000)) {
  444.                 k++;
  445.                 if (!(x & 0x40000000))
  446.                         return 32;
  447.                 }
  448.         return k;
  449.         }
  450.  
  451.  static int
  452. lo0bits
  453. #ifdef KR_headers
  454.         (y) unsigned long *y;
  455. #else
  456.         (unsigned long *y)
  457. #endif
  458. {
  459.         register int k;
  460.         register unsigned long x = *y;
  461.  
  462.         if (x & 7) {
  463.                 if (x & 1)
  464.                         return 0;
  465.                 if (x & 2) {
  466.                         *y = x >> 1;
  467.                         return 1;
  468.                         }
  469.                 *y = x >> 2;
  470.                 return 2;
  471.                 }
  472.         k = 0;
  473.         if (!(x & 0xffff)) {
  474.                 k = 16;
  475.                 x >>= 16;
  476.                 }
  477.         if (!(x & 0xff)) {
  478.                 k += 8;
  479.                 x >>= 8;
  480.                 }
  481.         if (!(x & 0xf)) {
  482.                 k += 4;
  483.                 x >>= 4;
  484.                 }
  485.         if (!(x & 0x3)) {
  486.                 k += 2;
  487.                 x >>= 2;
  488.                 }
  489.         if (!(x & 1)) {
  490.                 k++;
  491.                 x >>= 1;
  492.                 if (!x & 1)
  493.                         return 32;
  494.                 }
  495.         *y = x;
  496.         return k;
  497.         }
  498.  
  499.  static Bigint *
  500. i2b
  501. #ifdef KR_headers
  502.         (i) int i;
  503. #else
  504.         (int i)
  505. #endif
  506. {
  507.         Bigint *b;
  508.  
  509.         b = Balloc(1);
  510.         b->x[0] = i;
  511.         b->wds = 1;
  512.         return b;
  513.         }
  514.  
  515.  static Bigint *
  516. mult
  517. #ifdef KR_headers
  518.         (a, b) Bigint *a, *b;
  519. #else
  520.         (Bigint *a, Bigint *b)
  521. #endif
  522. {
  523.         Bigint *c;
  524.         int k, wa, wb, wc;
  525.         unsigned long carry, y, z;
  526.         unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  527. #ifdef Pack_32
  528.         unsigned long z2;
  529. #endif
  530.  
  531.         if (a->wds < b->wds) {
  532.                 c = a;
  533.                 a = b;
  534.                 b = c;
  535.                 }
  536.         k = a->k;
  537.         wa = a->wds;
  538.         wb = b->wds;
  539.         wc = wa + wb;
  540.         if (wc > a->maxwds)
  541.                 k++;
  542.         c = Balloc(k);
  543.         for(x = c->x, xa = x + wc; x < xa; x++)
  544.                 *x = 0;
  545.         xa = a->x;
  546.         xae = xa + wa;
  547.         xb = b->x;
  548.         xbe = xb + wb;
  549.         xc0 = c->x;
  550. #ifdef Pack_32
  551.         for(; xb < xbe; xb++, xc0++) {
  552.                 if (y = *xb & 0xffff) {
  553.                         x = xa;
  554.                         xc = xc0;
  555.                         carry = 0;
  556.                         do {
  557.                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  558.                                 carry = z >> 16;
  559.                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  560.                                 carry = z2 >> 16;
  561.                                 Storeinc(xc, z2, z);
  562.                                 }
  563.                                 while(x < xae);
  564.                         *xc = carry;
  565.                         }
  566.                 if (y = *xb >> 16) {
  567.                         x = xa;
  568.                         xc = xc0;
  569.                         carry = 0;
  570.                         z2 = *xc;
  571.                         do {
  572.                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  573.                                 carry = z >> 16;
  574.                                 Storeinc(xc, z, z2);
  575.                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  576.                                 carry = z2 >> 16;
  577.                                 }
  578.                                 while(x < xae);
  579.                         *xc = z2;
  580.                         }
  581.                 }
  582. #else
  583.         for(; xb < xbe; xc0++) {
  584.                 if (y = *xb++) {
  585.                         x = xa;
  586.                         xc = xc0;
  587.                         carry = 0;
  588.                         do {
  589.                                 z = *x++ * y + *xc + carry;
  590.                                 carry = z >> 16;
  591.                                 *xc++ = z & 0xffff;
  592.                                 }
  593.                                 while(x < xae);
  594.                         *xc = carry;
  595.                         }
  596.                 }
  597. #endif
  598.         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  599.         c->wds = wc;
  600.         return c;
  601.         }
  602.  
  603.  static Bigint *p5s;
  604.  
  605.  static Bigint *
  606. pow5mult
  607. #ifdef KR_headers
  608.         (b, k) Bigint *b; int k;
  609. #else
  610.         (Bigint *b, int k)
  611. #endif
  612. {
  613.         Bigint *b1, *p5, *p51;
  614.         int i;
  615.         static int p05[3] = { 5, 25, 125 };
  616.  
  617.         if (i = k & 3)
  618.                 b = multadd(b, p05[i-1], 0);
  619.  
  620.         if (!(k >>= 2))
  621.                 return b;
  622.         if (!(p5 = p5s)) {
  623.                 /* first time */
  624.                 p5 = p5s = i2b(625);
  625.                 p5->next = 0;
  626.                 }
  627.         for(;;) {
  628.                 if (k & 1) {
  629.                         b1 = mult(b, p5);
  630.                         Bfree(b);
  631.                         b = b1;
  632.                         }
  633.                 if (!(k >>= 1))
  634.                         break;
  635.                 if (!(p51 = p5->next)) {
  636.                         p51 = p5->next = mult(p5,p5);
  637.                         p51->next = 0;
  638.                         }
  639.                 p5 = p51;
  640.                 }
  641.         return b;
  642.         }
  643.  
  644.  static Bigint *
  645. lshift
  646. #ifdef KR_headers
  647.         (b, k) Bigint *b; int k;
  648. #else
  649.         (Bigint *b, int k)
  650. #endif
  651. {
  652.         int i, k1, n, n1;
  653.         Bigint *b1;
  654.         unsigned long *x, *x1, *xe, z;
  655.  
  656. #ifdef Pack_32
  657.         n = k >> 5;
  658. #else
  659.         n = k >> 4;
  660. #endif
  661.         k1 = b->k;
  662.         n1 = n + b->wds + 1;
  663.         for(i = b->maxwds; n1 > i; i <<= 1)
  664.                 k1++;
  665.         b1 = Balloc(k1);
  666.         x1 = b1->x;
  667.         for(i = 0; i < n; i++)
  668.                 *x1++ = 0;
  669.         x = b->x;
  670.         xe = x + b->wds;
  671. #ifdef Pack_32
  672.         if (k &= 0x1f) {
  673.                 k1 = 32 - k;
  674.                 z = 0;
  675.                 do {
  676.                         *x1++ = *x << k | z;
  677.                         z = *x++ >> k1;
  678.                         }
  679.                         while(x < xe);
  680.                 if (*x1 = z)
  681.                         ++n1;
  682.                 }
  683. #else
  684.         if (k &= 0xf) {
  685.                 k1 = 16 - k;
  686.                 z = 0;
  687.                 do {
  688.                         *x1++ = *x << k  & 0xffff | z;
  689.                         z = *x++ >> k1;
  690.                         }
  691.                         while(x < xe);
  692.                 if (*x1 = z)
  693.                         ++n1;
  694.                 }
  695. #endif
  696.         else do
  697.                 *x1++ = *x++;
  698.                 while(x < xe);
  699.         b1->wds = n1 - 1;
  700.         Bfree(b);
  701.         return b1;
  702.         }
  703.  
  704.  static int
  705. cmp
  706. #ifdef KR_headers
  707.         (a, b) Bigint *a, *b;
  708. #else
  709.         (Bigint *a, Bigint *b)
  710. #endif
  711. {
  712.         unsigned long *xa, *xa0, *xb, *xb0;
  713.         int i, j;
  714.  
  715.         i = a->wds;
  716.         j = b->wds;
  717. #ifdef DEBUG
  718.         if (i > 1 && !a->x[i-1])
  719.                 Bug("cmp called with a->x[a->wds-1] == 0");
  720.         if (j > 1 && !b->x[j-1])
  721.                 Bug("cmp called with b->x[b->wds-1] == 0");
  722. #endif
  723.         if (i -= j)
  724.                 return i;
  725.         xa0 = a->x;
  726.         xa = xa0 + j;
  727.         xb0 = b->x;
  728.         xb = xb0 + j;
  729.         for(;;) {
  730.                 if (*--xa != *--xb)
  731.                         return *xa < *xb ? -1 : 1;
  732.                 if (xa <= xa0)
  733.                         break;
  734.                 }
  735.         return 0;
  736.         }
  737.  
  738.  static Bigint *
  739. diff
  740. #ifdef KR_headers
  741.         (a, b) Bigint *a, *b;
  742. #else
  743.         (Bigint *a, Bigint *b)
  744. #endif
  745. {
  746.         Bigint *c;
  747.         int i, wa, wb;
  748.         long borrow, y; /* We need signed shifts here. */
  749.         unsigned long *xa, *xae, *xb, *xbe, *xc;
  750. #ifdef Pack_32
  751.         long z;
  752. #endif
  753.  
  754.         i = cmp(a,b);
  755.         if (!i) {
  756.                 c = Balloc(0);
  757.                 c->wds = 1;
  758.                 c->x[0] = 0;
  759.                 return c;
  760.                 }
  761.         if (i < 0) {
  762.                 c = a;
  763.                 a = b;
  764.                 b = c;
  765.                 i = 1;
  766.                 }
  767.         else
  768.                 i = 0;
  769.         c = Balloc(a->k);
  770.         c->sign = i;
  771.         wa = a->wds;
  772.         xa = a->x;
  773.         xae = xa + wa;
  774.         wb = b->wds;
  775.         xb = b->x;
  776.         xbe = xb + wb;
  777.         xc = c->x;
  778.         borrow = 0;
  779. #ifdef Pack_32
  780.         do {
  781.                 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
  782.                 borrow = y >> 16;
  783.                 Sign_Extend(borrow, y);
  784.                 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
  785.                 borrow = z >> 16;
  786.                 Sign_Extend(borrow, z);
  787.                 Storeinc(xc, z, y);
  788.                 }
  789.                 while(xb < xbe);
  790.         while(xa < xae) {
  791.                 y = (*xa & 0xffff) + borrow;
  792.                 borrow = y >> 16;
  793.                 Sign_Extend(borrow, y);
  794.                 z = (*xa++ >> 16) + borrow;
  795.                 borrow = z >> 16;
  796.                 Sign_Extend(borrow, z);
  797.                 Storeinc(xc, z, y);
  798.                 }
  799. #else
  800.         do {
  801.                 y = *xa++ - *xb++ + borrow;
  802.                 borrow = y >> 16;
  803.                 Sign_Extend(borrow, y);
  804.                 *xc++ = y & 0xffff;
  805.                 }
  806.                 while(xb < xbe);
  807.         while(xa < xae) {
  808.                 y = *xa++ + borrow;
  809.                 borrow = y >> 16;
  810.                 Sign_Extend(borrow, y);
  811.                 *xc++ = y & 0xffff;
  812.                 }
  813. #endif
  814.         while(!*--xc)
  815.                 wa--;
  816.         c->wds = wa;
  817.         return c;
  818.         }
  819.  
  820.  static double
  821. ulp
  822. #ifdef KR_headers
  823.         (x) double x;
  824. #else
  825.         (double x)
  826. #endif
  827. {
  828.         register long L;
  829.         double a;
  830.  
  831.         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
  832. #ifndef Sudden_Underflow
  833.         if (L > 0) {
  834. #endif
  835. #ifdef IBM
  836.                 L |= Exp_msk1 >> 4;
  837. #endif
  838.                 word0(a) = L;
  839.                 word1(a) = 0;
  840. #ifndef Sudden_Underflow
  841.                 }
  842.         else {
  843.                 L = -L >> Exp_shift;
  844.                 if (L < Exp_shift) {
  845.                         word0(a) = 0x80000 >> L;
  846.                         word1(a) = 0;
  847.                         }
  848.                 else {
  849.                         word0(a) = 0;
  850.                         L -= Exp_shift;
  851.                         word1(a) = L >= 31 ? 1 : 1 << 31 - L;
  852.                         }
  853.                 }
  854. #endif
  855.         return a;
  856.         }
  857.  
  858.  static double
  859. b2d
  860. #ifdef KR_headers
  861.         (a, e) Bigint *a; int *e;
  862. #else
  863.         (Bigint *a, int *e)
  864. #endif
  865. {
  866.         unsigned long *xa, *xa0, w, y, z;
  867.         int k;
  868.         double d;
  869. #ifdef VAX
  870.         unsigned long d0, d1;
  871. #else
  872. #define d0 word0(d)
  873. #define d1 word1(d)
  874. #endif
  875.  
  876.         xa0 = a->x;
  877.         xa = xa0 + a->wds;
  878.         y = *--xa;
  879. #ifdef DEBUG
  880.         if (!y) Bug("zero y in b2d");
  881. #endif
  882.         k = hi0bits(y);
  883.         *e = 32 - k;
  884. #ifdef Pack_32
  885.         if (k < Ebits) {
  886.                 d0 = Exp_1 | y >> Ebits - k;
  887.                 w = xa > xa0 ? *--xa : 0;
  888.                 d1 = y << (32-Ebits) + k | w >> Ebits - k;
  889.                 goto ret_d;
  890.                 }
  891.         z = xa > xa0 ? *--xa : 0;
  892.         if (k -= Ebits) {
  893.                 d0 = Exp_1 | y << k | z >> 32 - k;
  894.                 y = xa > xa0 ? *--xa : 0;
  895.                 d1 = z << k | y >> 32 - k;
  896.                 }
  897.         else {
  898.                 d0 = Exp_1 | y;
  899.                 d1 = z;
  900.                 }
  901. #else
  902.         if (k < Ebits + 16) {
  903.                 z = xa > xa0 ? *--xa : 0;
  904.                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  905.                 w = xa > xa0 ? *--xa : 0;
  906.                 y = xa > xa0 ? *--xa : 0;
  907.                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  908.                 goto ret_d;
  909.                 }
  910.         z = xa > xa0 ? *--xa : 0;
  911.         w = xa > xa0 ? *--xa : 0;
  912.         k -= Ebits + 16;
  913.         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  914.         y = xa > xa0 ? *--xa : 0;
  915.         d1 = w << k + 16 | y << k;
  916. #endif
  917.  ret_d:
  918. #ifdef VAX
  919.         word0(d) = d0 >> 16 | d0 << 16;
  920.         word1(d) = d1 >> 16 | d1 << 16;
  921. #else
  922. #undef d0
  923. #undef d1
  924. #endif
  925.         return d;
  926.         }
  927.  
  928.  static Bigint *
  929. d2b
  930. #ifdef KR_headers
  931.         (d, e, bits) double d; int *e, *bits;
  932. #else
  933.         (double d, int *e, int *bits)
  934. #endif
  935. {
  936.         Bigint *b;
  937.         int de, i, k;
  938.         unsigned long *x, y, z;
  939. #ifdef VAX
  940.         unsigned long d0, d1;
  941.         d0 = word0(d) >> 16 | word0(d) << 16;
  942.         d1 = word1(d) >> 16 | word1(d) << 16;
  943. #else
  944. #define d0 word0(d)
  945. #define d1 word1(d)
  946. #endif
  947.  
  948. #ifdef Pack_32
  949.         b = Balloc(1);
  950. #else
  951.         b = Balloc(2);
  952. #endif
  953.         x = b->x;
  954.  
  955.         z = d0 & Frac_mask;
  956.         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
  957. #ifdef Sudden_Underflow
  958.         de = (int)(d0 >> Exp_shift);
  959. #ifndef IBM
  960.         z |= Exp_msk11;
  961. #endif
  962. #else
  963.         if (de = (int)(d0 >> Exp_shift))
  964.                 z |= Exp_msk1;
  965. #endif
  966. #ifdef Pack_32
  967.         if (y = d1) {
  968.                 if (k = lo0bits(&y)) {
  969.                         x[0] = y | z << 32 - k;
  970.                         z >>= k;
  971.                         }
  972.                 else
  973.                         x[0] = y;
  974.                 i = b->wds = (x[1] = z) ? 2 : 1;
  975.                 }
  976.         else {
  977. #ifdef DEBUG
  978.                 if (!z)
  979.                         Bug("Zero passed to d2b");
  980. #endif
  981.                 k = lo0bits(&z);
  982.                 x[0] = z;
  983.                 i = b->wds = 1;
  984.                 k += 32;
  985.                 }
  986. #else
  987.         if (y = d1) {
  988.                 if (k = lo0bits(&y))
  989.                         if (k >= 16) {
  990.                                 x[0] = y | z << 32 - k & 0xffff;
  991.                                 x[1] = z >> k - 16 & 0xffff;
  992.                                 x[2] = z >> k;
  993.                                 i = 2;
  994.                                 }
  995.                         else {
  996.                                 x[0] = y & 0xffff;
  997.                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
  998.                                 x[2] = z >> k & 0xffff;
  999.                                 x[3] = z >> k+16;
  1000.                                 i = 3;
  1001.                                 }
  1002.                 else {
  1003.                         x[0] = y & 0xffff;
  1004.                         x[1] = y >> 16;
  1005.                         x[2] = z & 0xffff;
  1006.                         x[3] = z >> 16;
  1007.                         i = 3;
  1008.                         }
  1009.                 }
  1010.         else {
  1011. #ifdef DEBUG
  1012.                 if (!z)
  1013.                         Bug("Zero passed to d2b");
  1014. #endif
  1015.                 k = lo0bits(&z);
  1016.                 if (k >= 16) {
  1017.                         x[0] = z;
  1018.                         i = 0;
  1019.                         }
  1020.                 else {
  1021.                         x[0] = z & 0xffff;
  1022.                         x[1] = z >> 16;
  1023.                         i = 1;
  1024.                         }
  1025.                 k += 32;
  1026.                 }
  1027.         while(!x[i])
  1028.                 --i;
  1029.         b->wds = i + 1;
  1030. #endif
  1031. #ifndef Sudden_Underflow
  1032.         if (de) {
  1033. #endif
  1034. #ifdef IBM
  1035.                 *e = (de - Bias - (P-1) << 2) + k;
  1036.                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
  1037. #else
  1038.                 *e = de - Bias - (P-1) + k;
  1039.                 *bits = P - k;
  1040. #endif
  1041. #ifndef Sudden_Underflow
  1042.                 }
  1043.         else {
  1044.                 *e = de - Bias - (P-1) + 1 + k;
  1045. #ifdef Pack_32
  1046.                 *bits = 32*i - hi0bits(x[i-1]);
  1047. #else
  1048.                 *bits = (i+2)*16 - hi0bits(x[i]);
  1049. #endif
  1050.                 }
  1051. #endif
  1052.         return b;
  1053.         }
  1054. #undef d0
  1055. #undef d1
  1056.  
  1057.  static double
  1058. ratio
  1059. #ifdef KR_headers
  1060.         (a, b) Bigint *a, *b;
  1061. #else
  1062.         (Bigint *a, Bigint *b)
  1063. #endif
  1064. {
  1065.         double da, db;
  1066.         int k, ka, kb;
  1067.  
  1068.         da = b2d(a, &ka);
  1069.         db = b2d(b, &kb);
  1070. #ifdef Pack_32
  1071.         k = ka - kb + 32*(a->wds - b->wds);
  1072. #else
  1073.         k = ka - kb + 16*(a->wds - b->wds);
  1074. #endif
  1075. #ifdef IBM
  1076.         if (k > 0) {
  1077.                 word0(da) += (k >> 2)*Exp_msk1;
  1078.                 if (k &= 3)
  1079.                         da *= 1 << k;
  1080.                 }
  1081.         else {
  1082.                 k = -k;
  1083.                 word0(db) += (k >> 2)*Exp_msk1;
  1084.                 if (k &= 3)
  1085.                         db *= 1 << k;
  1086.                 }
  1087. #else
  1088.         if (k > 0)
  1089.                 word0(da) += k*Exp_msk1;
  1090.         else {
  1091.                 k = -k;
  1092.                 word0(db) += k*Exp_msk1;
  1093.                 }
  1094. #endif
  1095.         return da / db;
  1096.         }
  1097.  
  1098.  static double
  1099. tens[] = {
  1100.                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  1101.                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  1102.                 1e20, 1e21, 1e22
  1103. #ifdef VAX
  1104.                 , 1e23, 1e24
  1105. #endif
  1106.                 };
  1107.  
  1108.  static double
  1109. #ifdef IEEE_Arith
  1110. bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  1111. static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
  1112. #define n_bigtens 5
  1113. #else
  1114. #ifdef IBM
  1115. bigtens[] = { 1e16, 1e32, 1e64 };
  1116. static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  1117. #define n_bigtens 3
  1118. #else
  1119. bigtens[] = { 1e16, 1e32 };
  1120. static double tinytens[] = { 1e-16, 1e-32 };
  1121. #define n_bigtens 2
  1122. #endif
  1123. #endif
  1124.  
  1125.  double
  1126. strtod
  1127. #ifdef KR_headers
  1128.         (s00, se) CONST char *s00; char **se;
  1129. #else
  1130.         (CONST char *s00, char **se)
  1131. #endif
  1132. {
  1133.         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
  1134.                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
  1135.         CONST char *s, *s0, *s1;
  1136.         double aadj, aadj1, adj, rv, rv0;
  1137.         long L;
  1138.         unsigned long y, z;
  1139.         Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
  1140.     TEST_ENDIANNESS;
  1141.         sign = nz0 = nz = 0;
  1142.         rv = 0.;
  1143.         for(s = s00;;s++) switch(*s) {
  1144.                 case '-':
  1145.                         sign = 1;
  1146.                         /* no break */
  1147.                 case '+':
  1148.                         if (*++s)
  1149.                                 goto break2;
  1150.                         /* no break */
  1151.                 case 0:
  1152.                         goto ret;
  1153.                 case '\t':
  1154.                 case '\n':
  1155.                 case '\v':
  1156.                 case '\f':
  1157.                 case '\r':
  1158.                 case ' ':
  1159.                         continue;
  1160.                 default:
  1161.                         goto break2;
  1162.                 }
  1163.  break2:
  1164.         if (*s == '0') {
  1165.                 nz0 = 1;
  1166.                 while(*++s == '0') ;
  1167.                 if (!*s)
  1168.                         goto ret;
  1169.                 }
  1170.         s0 = s;
  1171.         y = z = 0;
  1172.         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  1173.                 if (nd < 9)
  1174.                         y = 10*y + c - '0';
  1175.                 else if (nd < 16)
  1176.                         z = 10*z + c - '0';
  1177.         nd0 = nd;
  1178.         if (c == '.') {
  1179.                 c = *++s;
  1180.                 if (!nd) {
  1181.                         for(; c == '0'; c = *++s)
  1182.                                 nz++;
  1183.                         if (c > '0' && c <= '9') {
  1184.                                 s0 = s;
  1185.                                 nf += nz;
  1186.                                 nz = 0;
  1187.                                 goto have_dig;
  1188.                                 }
  1189.                         goto dig_done;
  1190.                         }
  1191.                 for(; c >= '0' && c <= '9'; c = *++s) {
  1192.  have_dig:
  1193.                         nz++;
  1194.                         if (c -= '0') {
  1195.                                 nf += nz;
  1196.                                 for(i = 1; i < nz; i++)
  1197.                                         if (nd++ < 9)
  1198.                                                 y *= 10;
  1199.                                         else if (nd <= DBL_DIG + 1)
  1200.                                                 z *= 10;
  1201.                                 if (nd++ < 9)
  1202.                                         y = 10*y + c;
  1203.                                 else if (nd <= DBL_DIG + 1)
  1204.                                         z = 10*z + c;
  1205.                                 nz = 0;
  1206.                                 }
  1207.                         }
  1208.                 }
  1209.  dig_done:
  1210.         e = 0;
  1211.         if (c == 'e' || c == 'E') {
  1212.                 if (!nd && !nz && !nz0) {
  1213.                         s = s00;
  1214.                         goto ret;
  1215.                         }
  1216.                 s00 = s;
  1217.                 esign = 0;
  1218.                 switch(c = *++s) {
  1219.                         case '-':
  1220.                                 esign = 1;
  1221.                         case '+':
  1222.                                 c = *++s;
  1223.                         }
  1224.                 if (c >= '0' && c <= '9') {
  1225.                         while(c == '0')
  1226.                                 c = *++s;
  1227.                         if (c > '0' && c <= '9') {
  1228.                                 e = c - '0';
  1229.                                 s1 = s;
  1230.                                 while((c = *++s) >= '0' && c <= '9')
  1231.                                         e = 10*e + c - '0';
  1232.                                 if (s - s1 > 8)
  1233.                                         /* Avoid confusion from exponents
  1234.                                          * so large that e might overflow.
  1235.                                          */
  1236.                                         e = 9999999;
  1237.                                 if (esign)
  1238.                                         e = -e;
  1239.                                 }
  1240.                         else
  1241.                                 e = 0;
  1242.                         }
  1243.                 else
  1244.                         s = s00;
  1245.                 }
  1246.         if (!nd) {
  1247.                 if (!nz && !nz0)
  1248.                         s = s00;
  1249.                 goto ret;
  1250.                 }
  1251.         e1 = e -= nf;
  1252.  
  1253.         /* Now we have nd0 digits, starting at s0, followed by a
  1254.          * decimal point, followed by nd-nd0 digits.  The number we're
  1255.          * after is the integer represented by those digits times
  1256.          * 10**e */
  1257.  
  1258.         if (!nd0)
  1259.                 nd0 = nd;
  1260.         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  1261.         rv = y;
  1262.         if (k > 9)
  1263.                 rv = tens[k - 9] * rv + z;
  1264.         if (nd <= DBL_DIG
  1265. #ifndef RND_PRODQUOT
  1266.                 && FLT_ROUNDS == 1
  1267. #endif
  1268.                         ) {
  1269.                 if (!e)
  1270.                         goto ret;
  1271.                 if (e > 0) {
  1272.                         if (e <= Ten_pmax) {
  1273. #ifdef VAX
  1274.                                 goto vax_ovfl_check;
  1275. #else
  1276.                                 /* rv = */ rounded_product(rv, tens[e]);
  1277.                                 goto ret;
  1278. #endif
  1279.                                 }
  1280.                         i = DBL_DIG - nd;
  1281.                         if (e <= Ten_pmax + i) {
  1282.                                 /* A fancier test would sometimes let us do
  1283.                                  * this for larger i values.
  1284.                                  */
  1285.                                 e -= i;
  1286.                                 rv *= tens[i];
  1287. #ifdef VAX
  1288.                                 /* VAX exponent range is so narrow we must
  1289.                                  * worry about overflow here...
  1290.                                  */
  1291.  vax_ovfl_check:
  1292.                                 word0(rv) -= P*Exp_msk1;
  1293.                                 /* rv = */ rounded_product(rv, tens[e]);
  1294.                                 if ((word0(rv) & Exp_mask)
  1295.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  1296.                                         goto ovfl;
  1297.                                 word0(rv) += P*Exp_msk1;
  1298. #else
  1299.                                 /* rv = */ rounded_product(rv, tens[e]);
  1300. #endif
  1301.                                 goto ret;
  1302.                                 }
  1303.                         }
  1304. #ifndef Inaccurate_Divide
  1305.                 else if (e >= -Ten_pmax) {
  1306.                         /* rv = */ rounded_quotient(rv, tens[-e]);
  1307.                         goto ret;
  1308.                         }
  1309. #endif
  1310.                 }
  1311.         e1 += nd - k;
  1312.  
  1313.         /* Get starting approximation = rv * 10**e1 */
  1314.  
  1315.         if (e1 > 0) {
  1316.                 if (i = e1 & 15)
  1317.                         rv *= tens[i];
  1318.                 if (e1 &= ~15) {
  1319.                         if (e1 > DBL_MAX_10_EXP) {
  1320.  ovfl:
  1321.                                 errno = ERANGE;
  1322. #ifndef HUGE_VAL
  1323. #define HUGE_VAL        1.7976931348623157E+308
  1324. #endif
  1325.                                 rv = HUGE_VAL;
  1326.                                 goto ret;
  1327.                                 }
  1328.                         if (e1 >>= 4) {
  1329.                                 for(j = 0; e1 > 1; j++, e1 >>= 1)
  1330.                                         if (e1 & 1)
  1331.                                                 rv *= bigtens[j];
  1332.                         /* The last multiplication could overflow. */
  1333.                                 word0(rv) -= P*Exp_msk1;
  1334.                                 rv *= bigtens[j];
  1335.                                 if ((z = word0(rv) & Exp_mask)
  1336.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  1337.                                         goto ovfl;
  1338.                                 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  1339.                                         /* set to largest number */
  1340.                                         /* (Can't trust DBL_MAX) */
  1341.                                         word0(rv) = Big0;
  1342.                                         word1(rv) = Big1;
  1343.                                         }
  1344.                                 else
  1345.                                         word0(rv) += P*Exp_msk1;
  1346.                                 }
  1347.  
  1348.                         }
  1349.                 }
  1350.         else if (e1 < 0) {
  1351.                 e1 = -e1;
  1352.                 if (i = e1 & 15)
  1353.                         rv /= tens[i];
  1354.                 if (e1 &= ~15) {
  1355.                         e1 >>= 4;
  1356.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  1357.                                 if (e1 & 1)
  1358.                                         rv *= tinytens[j];
  1359.                         /* The last multiplication could underflow. */
  1360.                         rv0 = rv;
  1361.                         rv *= tinytens[j];
  1362.                         if (!rv) {
  1363.                                 rv = 2.*rv0;
  1364.                                 rv *= tinytens[j];
  1365.                                 if (!rv) {
  1366.  undfl:
  1367.                                         rv = 0.;
  1368.                                         errno = ERANGE;
  1369.                                         goto ret;
  1370.                                         }
  1371.                                 word0(rv) = Tiny0;
  1372.                                 word1(rv) = Tiny1;
  1373.                                 /* The refinement below will clean
  1374.                                  * this approximation up.
  1375.                                  */
  1376.                                 }
  1377.                         }
  1378.                 }
  1379.  
  1380.         /* Now the hard part -- adjusting rv to the correct value.*/
  1381.  
  1382.         /* Put digits into bd: true value = bd * 10^e */
  1383.  
  1384.         bd0 = s2b(s0, nd0, nd, y);
  1385.  
  1386.         for(;;) {
  1387.                 bd = Balloc(bd0->k);
  1388.                 Bcopy(bd, bd0);
  1389.                 bb = d2b(rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
  1390.                 bs = i2b(1);
  1391.  
  1392.                 if (e >= 0) {
  1393.                         bb2 = bb5 = 0;
  1394.                         bd2 = bd5 = e;
  1395.                         }
  1396.                 else {
  1397.                         bb2 = bb5 = -e;
  1398.                         bd2 = bd5 = 0;
  1399.                         }
  1400.                 if (bbe >= 0)
  1401.                         bb2 += bbe;
  1402.                 else
  1403.                         bd2 -= bbe;
  1404.                 bs2 = bb2;
  1405. #ifdef Sudden_Underflow
  1406. #ifdef IBM
  1407.                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  1408. #else
  1409.                 j = P + 1 - bbbits;
  1410. #endif
  1411. #else
  1412.                 i = bbe + bbbits - 1;   /* logb(rv) */
  1413.                 if (i < Emin)   /* denormal */
  1414.                         j = bbe + (P-Emin);
  1415.                 else
  1416.                         j = P + 1 - bbbits;
  1417. #endif
  1418.                 bb2 += j;
  1419.                 bd2 += j;
  1420.                 i = bb2 < bd2 ? bb2 : bd2;
  1421.                 if (i > bs2)
  1422.                         i = bs2;
  1423.                 if (i > 0) {
  1424.                         bb2 -= i;
  1425.                         bd2 -= i;
  1426.                         bs2 -= i;
  1427.                         }
  1428.                 if (bb5 > 0) {
  1429.                         bs = pow5mult(bs, bb5);
  1430.                         bb1 = mult(bs, bb);
  1431.                         Bfree(bb);
  1432.                         bb = bb1;
  1433.                         }
  1434.                 if (bb2 > 0)
  1435.                         bb = lshift(bb, bb2);
  1436.                 if (bd5 > 0)
  1437.                         bd = pow5mult(bd, bd5);
  1438.                 if (bd2 > 0)
  1439.                         bd = lshift(bd, bd2);
  1440.                 if (bs2 > 0)
  1441.                         bs = lshift(bs, bs2);
  1442.                 delta = diff(bb, bd);
  1443.                 dsign = delta->sign;
  1444.                 delta->sign = 0;
  1445.                 i = cmp(delta, bs);
  1446.                 if (i < 0) {
  1447.                         /* Error is less than half an ulp -- check for
  1448.                          * special case of mantissa a power of two.
  1449.                          */
  1450.                         if (dsign || word1(rv) || word0(rv) & Bndry_mask)
  1451.                                 break;
  1452.                         delta = lshift(delta,Log2P);
  1453.                         if (cmp(delta, bs) > 0)
  1454.                                 goto drop_down;
  1455.                         break;
  1456.                         }
  1457.                 if (i == 0) {
  1458.                         /* exactly half-way between */
  1459.                         if (dsign) {
  1460.                                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
  1461.                                  &&  word1(rv) == 0xffffffff) {
  1462.                                         /*boundary case -- increment exponent*/
  1463.                                         word0(rv) = (word0(rv) & Exp_mask)
  1464.                                                 + Exp_msk1
  1465. #ifdef IBM
  1466.                                                 | Exp_msk1 >> 4
  1467. #endif
  1468.                                                 ;
  1469.                                         word1(rv) = 0;
  1470.                                         break;
  1471.                                         }
  1472.                                 }
  1473.                         else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
  1474.  drop_down:
  1475.                                 /* boundary case -- decrement exponent */
  1476. #ifdef Sudden_Underflow
  1477.                                 L = word0(rv) & Exp_mask;
  1478. #ifdef IBM
  1479.                                 if (L <  Exp_msk1)
  1480. #else
  1481.                                 if (L <= Exp_msk1)
  1482. #endif
  1483.                                         goto undfl;
  1484.                                 L -= Exp_msk1;
  1485. #else
  1486.                                 L = (word0(rv) & Exp_mask) - Exp_msk1;
  1487. #endif
  1488.                                 word0(rv) = L | Bndry_mask1;
  1489.                                 word1(rv) = 0xffffffff;
  1490. #ifdef IBM
  1491.                                 goto cont;
  1492. #else
  1493.                                 break;
  1494. #endif
  1495.                                 }
  1496. #ifndef ROUND_BIASED
  1497.                         if (!(word1(rv) & LSB))
  1498.                                 break;
  1499. #endif
  1500.                         if (dsign)
  1501.                                 rv += ulp(rv);
  1502. #ifndef ROUND_BIASED
  1503.                         else {
  1504.                                 rv -= ulp(rv);
  1505. #ifndef Sudden_Underflow
  1506.                                 if (!rv)
  1507.                                         goto undfl;
  1508. #endif
  1509.                                 }
  1510. #endif
  1511.                         break;
  1512.                         }
  1513.                 if ((aadj = ratio(delta, bs)) <= 2.) {
  1514.                         if (dsign)
  1515.                                 aadj = aadj1 = 1.;
  1516.                         else if (word1(rv) || word0(rv) & Bndry_mask) {
  1517. #ifndef Sudden_Underflow
  1518.                                 if (word1(rv) == Tiny1 && !word0(rv))
  1519.                                         goto undfl;
  1520. #endif
  1521.                                 aadj = 1.;
  1522.                                 aadj1 = -1.;
  1523.                                 }
  1524.                         else {
  1525.                                 /* special case -- power of FLT_RADIX to be */
  1526.                                 /* rounded down... */
  1527.  
  1528.                                 if (aadj < 2./FLT_RADIX)
  1529.                                         aadj = 1./FLT_RADIX;
  1530.                                 else
  1531.                                         aadj *= 0.5;
  1532.                                 aadj1 = -aadj;
  1533.                                 }
  1534.                         }
  1535.                 else {
  1536.                         aadj *= 0.5;
  1537.                         aadj1 = dsign ? aadj : -aadj;
  1538. #ifdef Check_FLT_ROUNDS
  1539.                         switch(FLT_ROUNDS) {
  1540.                                 case 2: /* towards +infinity */
  1541.                                         aadj1 -= 0.5;
  1542.                                         break;
  1543.                                 case 0: /* towards 0 */
  1544.                                 case 3: /* towards -infinity */
  1545.                                         aadj1 += 0.5;
  1546.                                 }
  1547. #else
  1548.                         if (FLT_ROUNDS == 0)
  1549.                                 aadj1 += 0.5;
  1550. #endif
  1551.                         }
  1552.                 y = word0(rv) & Exp_mask;
  1553.  
  1554.                 /* Check for overflow */
  1555.  
  1556.                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  1557.                         rv0 = rv;
  1558.                         word0(rv) -= P*Exp_msk1;
  1559.                         adj = aadj1 * ulp(rv);
  1560.                         rv += adj;
  1561.                         if ((word0(rv) & Exp_mask) >=
  1562.                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  1563.                                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
  1564.                                         goto ovfl;
  1565.                                 word0(rv) = Big0;
  1566.                                 word1(rv) = Big1;
  1567.                                 goto cont;
  1568.                                 }
  1569.                         else
  1570.                                 word0(rv) += P*Exp_msk1;
  1571.                         }
  1572.                 else {
  1573. #ifdef Sudden_Underflow
  1574.                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
  1575.                                 rv0 = rv;
  1576.                                 word0(rv) += P*Exp_msk1;
  1577.                                 adj = aadj1 * ulp(rv);
  1578.                                 rv += adj;
  1579. #ifdef IBM
  1580.                                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
  1581. #else
  1582.                                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
  1583. #endif
  1584.                                         {
  1585.                                         if (word0(rv0) == Tiny0
  1586.                                          && word1(rv0) == Tiny1)
  1587.                                                 goto undfl;
  1588.                                         word0(rv) = Tiny0;
  1589.                                         word1(rv) = Tiny1;
  1590.                                         goto cont;
  1591.                                         }
  1592.                                 else
  1593.                                         word0(rv) -= P*Exp_msk1;
  1594.                                 }
  1595.                         else {
  1596.                                 adj = aadj1 * ulp(rv);
  1597.                                 rv += adj;
  1598.                                 }
  1599. #else
  1600.                         /* Compute adj so that the IEEE rounding rules will
  1601.                          * correctly round rv + adj in some half-way cases.
  1602.                          * If rv * ulp(rv) is denormalized (i.e.,
  1603.                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  1604.                          * trouble from bits lost to denormalization;
  1605.                          * example: 1.2e-307 .
  1606.                          */
  1607.                         if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
  1608.                                 aadj1 = (double)(int)(aadj + 0.5);
  1609.                                 if (!dsign)
  1610.                                         aadj1 = -aadj1;
  1611.                                 }
  1612.                         adj = aadj1 * ulp(rv);
  1613.                         rv += adj;
  1614. #endif
  1615.                         }
  1616.                 z = word0(rv) & Exp_mask;
  1617.                 if (y == z) {
  1618.                         /* Can we stop now? */
  1619.                         L = (long)aadj;
  1620.                         aadj -= L;
  1621.                         /* The tolerances below are conservative. */
  1622.                         if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
  1623.                                 if (aadj < .4999999 || aadj > .5000001)
  1624.                                         break;
  1625.                                 }
  1626.                         else if (aadj < .4999999/FLT_RADIX)
  1627.                                 break;
  1628.                         }
  1629.  cont:
  1630.                 Bfree(bb);
  1631.                 Bfree(bd);
  1632.                 Bfree(bs);
  1633.                 Bfree(delta);
  1634.                 }
  1635.         Bfree(bb);
  1636.         Bfree(bd);
  1637.         Bfree(bs);
  1638.         Bfree(bd0);
  1639.         Bfree(delta);
  1640.  ret:
  1641.         if (se)
  1642.                 *se = (char *)s;
  1643.         return sign ? -rv : rv;
  1644.         }
  1645.  
  1646.  static int
  1647. quorem
  1648. #ifdef KR_headers
  1649.         (b, S) Bigint *b, *S;
  1650. #else
  1651.         (Bigint *b, Bigint *S)
  1652. #endif
  1653. {
  1654.         int n;
  1655.         long borrow, y;
  1656.         unsigned long carry, q, ys;
  1657.         unsigned long *bx, *bxe, *sx, *sxe;
  1658. #ifdef Pack_32
  1659.         long z;
  1660.         unsigned long si, zs;
  1661. #endif
  1662.  
  1663.         n = S->wds;
  1664. #ifdef DEBUG
  1665.         /*debug*/ if (b->wds > n)
  1666.         /*debug*/       Bug("oversize b in quorem");
  1667. #endif
  1668.         if (b->wds < n)
  1669.                 return 0;
  1670.         sx = S->x;
  1671.         sxe = sx + --n;
  1672.         bx = b->x;
  1673.         bxe = bx + n;
  1674.         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
  1675. #ifdef DEBUG
  1676.         /*debug*/ if (q > 9)
  1677.         /*debug*/       Bug("oversized quotient in quorem");
  1678. #endif
  1679.         if (q) {
  1680.                 borrow = 0;
  1681.                 carry = 0;
  1682.                 do {
  1683. #ifdef Pack_32
  1684.                         si = *sx++;
  1685.                         ys = (si & 0xffff) * q + carry;
  1686.                         zs = (si >> 16) * q + (ys >> 16);
  1687.                         carry = zs >> 16;
  1688.                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1689.                         borrow = y >> 16;
  1690.                         Sign_Extend(borrow, y);
  1691.                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1692.                         borrow = z >> 16;
  1693.                         Sign_Extend(borrow, z);
  1694.                         Storeinc(bx, z, y);
  1695. #else
  1696.                         ys = *sx++ * q + carry;
  1697.                         carry = ys >> 16;
  1698.                         y = *bx - (ys & 0xffff) + borrow;
  1699.                         borrow = y >> 16;
  1700.                         Sign_Extend(borrow, y);
  1701.                         *bx++ = y & 0xffff;
  1702. #endif
  1703.                         }
  1704.                         while(sx <= sxe);
  1705.                 if (!*bxe) {
  1706.                         bx = b->x;
  1707.                         while(--bxe > bx && !*bxe)
  1708.                                 --n;
  1709.                         b->wds = n;
  1710.                         }
  1711.                 }
  1712.         if (cmp(b, S) >= 0) {
  1713.                 q++;
  1714.                 borrow = 0;
  1715.                 carry = 0;
  1716.                 bx = b->x;
  1717.                 sx = S->x;
  1718.                 do {
  1719. #ifdef Pack_32
  1720.                         si = *sx++;
  1721.                         ys = (si & 0xffff) + carry;
  1722.                         zs = (si >> 16) + (ys >> 16);
  1723.                         carry = zs >> 16;
  1724.                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1725.                         borrow = y >> 16;
  1726.                         Sign_Extend(borrow, y);
  1727.                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1728.                         borrow = z >> 16;
  1729.                         Sign_Extend(borrow, z);
  1730.                         Storeinc(bx, z, y);
  1731. #else
  1732.                         ys = *sx++ + carry;
  1733.                         carry = ys >> 16;
  1734.                         y = *bx - (ys & 0xffff) + borrow;
  1735.                         borrow = y >> 16;
  1736.                         Sign_Extend(borrow, y);
  1737.                         *bx++ = y & 0xffff;
  1738. #endif
  1739.                         }
  1740.                         while(sx <= sxe);
  1741.                 bx = b->x;
  1742.                 bxe = bx + n;
  1743.                 if (!*bxe) {
  1744.                         while(--bxe > bx && !*bxe)
  1745.                                 --n;
  1746.                         b->wds = n;
  1747.                         }
  1748.                 }
  1749.         return q;
  1750.         }
  1751.  
  1752. /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  1753.  *
  1754.  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  1755.  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
  1756.  *
  1757.  * Modifications:
  1758.  *      1. Rather than iterating, we use a simple numeric overestimate
  1759.  *         to determine k = floor(log10(d)).  We scale relevant
  1760.  *         quantities using O(log2(k)) rather than O(k) multiplications.
  1761.  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  1762.  *         try to generate digits strictly left to right.  Instead, we
  1763.  *         compute with fewer bits and propagate the carry if necessary
  1764.  *         when rounding the final digit up.  This is often faster.
  1765.  *      3. Under the assumption that input will be rounded nearest,
  1766.  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  1767.  *         That is, we allow equality in stopping tests when the
  1768.  *         round-nearest rule will give the same floating-point value
  1769.  *         as would satisfaction of the stopping test with strict
  1770.  *         inequality.
  1771.  *      4. We remove common factors of powers of 2 from relevant
  1772.  *         quantities.
  1773.  *      5. When converting floating-point integers less than 1e16,
  1774.  *         we use floating-point arithmetic rather than resorting
  1775.  *         to multiple-precision integers.
  1776.  *      6. When asked to produce fewer than 15 digits, we first try
  1777.  *         to get by with floating-point arithmetic; we resort to
  1778.  *         multiple-precision integer arithmetic only if we cannot
  1779.  *         guarantee that the floating-point calculation has given
  1780.  *         the correctly rounded result.  For k requested digits and
  1781.  *         "uniformly" distributed input, the probability is
  1782.  *         something like 10^(k-15) that we must resort to the long
  1783.  *         calculation.
  1784.  */
  1785.  
  1786.  char *
  1787. dtoa
  1788. #ifdef KR_headers
  1789.         (d, mode, ndigits, decpt, sign, rve)
  1790.         double d; int mode, ndigits, *decpt, *sign; char **rve;
  1791. #else
  1792.         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  1793. #endif
  1794. {
  1795.  /*     Arguments ndigits, decpt, sign are similar to those
  1796.         of ecvt and fcvt; trailing zeros are suppressed from
  1797.         the returned string.  If not null, *rve is set to point
  1798.         to the end of the return value.  If d is +-Infinity or NaN,
  1799.         then *decpt is set to 9999.
  1800.  
  1801.         mode:
  1802.                 0 ==> shortest string that yields d when read in
  1803.                         and rounded to nearest.
  1804.                 1 ==> like 0, but with Steele & White stopping rule;
  1805.                         e.g. with IEEE P754 arithmetic , mode 0 gives
  1806.                         1e23 whereas mode 1 gives 9.999999999999999e22.
  1807.                 2 ==> max(1,ndigits) significant digits.  This gives a
  1808.                         return value similar to that of ecvt, except
  1809.                         that trailing zeros are suppressed.
  1810.                 3 ==> through ndigits past the decimal point.  This
  1811.                         gives a return value similar to that from fcvt,
  1812.                         except that trailing zeros are suppressed, and
  1813.                         ndigits can be negative.
  1814.                 4-9 should give the same return values as 2-3, i.e.,
  1815.                         4 <= mode <= 9 ==> same return as mode
  1816.                         2 + (mode & 1).  These modes are mainly for
  1817.                         debugging; often they run slower but sometimes
  1818.                         faster than modes 2-3.
  1819.                 4,5,8,9 ==> left-to-right digit generation.
  1820.                 6-9 ==> don't try fast floating-point estimate
  1821.                         (if applicable).
  1822.  
  1823.                 Values of mode other than 0-9 are treated as mode 0.
  1824.  
  1825.                 Sufficient space is allocated to the return value
  1826.                 to hold the suppressed trailing zeros.
  1827.         */
  1828.  
  1829.         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
  1830.                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
  1831.                 spec_case, try_quick;
  1832.         long L;
  1833. #ifndef Sudden_Underflow
  1834.         int denorm;
  1835.         unsigned long x;
  1836. #endif
  1837.         Bigint *b, *b1, *delta, *mlo, *mhi, *S;
  1838.         double d2, ds, eps;
  1839.         char *s, *s0;
  1840.         static Bigint *result;
  1841.         static int result_k;
  1842.  
  1843.     TEST_ENDIANNESS;
  1844.         if (result) {
  1845.                 result->k = result_k;
  1846.                 result->maxwds = 1 << result_k;
  1847.                 Bfree(result);
  1848.                 result = 0;
  1849.                 }
  1850.  
  1851.         if (word0(d) & Sign_bit) {
  1852.                 /* set sign for everything, including 0's and NaNs */
  1853.                 *sign = 1;
  1854.                 word0(d) &= ~Sign_bit;  /* clear sign bit */
  1855.                 }
  1856.         else
  1857.                 *sign = 0;
  1858.  
  1859. #if defined(IEEE_Arith) + defined(VAX)
  1860. #ifdef IEEE_Arith
  1861.         if ((word0(d) & Exp_mask) == Exp_mask)
  1862. #else
  1863.         if (word0(d)  == 0x8000)
  1864. #endif
  1865.                 {
  1866.                 /* Infinity or NaN */
  1867.                 *decpt = 9999;
  1868.                 s =
  1869. #ifdef IEEE_Arith
  1870.                         !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
  1871. #endif
  1872.                                 "NaN";
  1873.                 if (rve)
  1874.                         *rve =
  1875. #ifdef IEEE_Arith
  1876.                                 s[3] ? s + 8 :
  1877. #endif
  1878.                                                 s + 3;
  1879.                 return s;
  1880.                 }
  1881. #endif
  1882. #ifdef IBM
  1883.         d += 0; /* normalize */
  1884. #endif
  1885.         if (!d) {
  1886.                 *decpt = 1;
  1887.                 s = "0";
  1888.                 if (rve)
  1889.                         *rve = s + 1;
  1890.                 return s;
  1891.                 }
  1892.  
  1893.         b = d2b(d, &be, &bbits);
  1894. #ifdef Sudden_Underflow
  1895.         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
  1896. #else
  1897.         if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
  1898. #endif
  1899.                 d2 = d;
  1900.                 word0(d2) &= Frac_mask1;
  1901.                 word0(d2) |= Exp_11;
  1902. #ifdef IBM
  1903.                 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
  1904.                         d2 /= 1 << j;
  1905. #endif
  1906.  
  1907.                 /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
  1908.                  * log10(x)      =  log(x) / log(10)
  1909.                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
  1910.                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
  1911.                  *
  1912.                  * This suggests computing an approximation k to log10(d) by
  1913.                  *
  1914.                  * k = (i - Bias)*0.301029995663981
  1915.                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
  1916.                  *
  1917.                  * We want k to be too large rather than too small.
  1918.                  * The error in the first-order Taylor series approximation
  1919.                  * is in our favor, so we just round up the constant enough
  1920.                  * to compensate for any error in the multiplication of
  1921.                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
  1922.                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
  1923.                  * adding 1e-13 to the constant term more than suffices.
  1924.                  * Hence we adjust the constant term to 0.1760912590558.
  1925.                  * (We could get a more accurate k by invoking log10,
  1926.                  *  but this is probably not worthwhile.)
  1927.                  */
  1928.  
  1929.                 i -= Bias;
  1930. #ifdef IBM
  1931.                 i <<= 2;
  1932.                 i += j;
  1933. #endif
  1934. #ifndef Sudden_Underflow
  1935.                 denorm = 0;
  1936.                 }
  1937.         else {
  1938.                 /* d is denormalized */
  1939.  
  1940.                 i = bbits + be + (Bias + (P-1) - 1);
  1941.                 x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
  1942.                             : word1(d) << 32 - i;
  1943.                 d2 = x;
  1944.                 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
  1945.                 i -= (Bias + (P-1) - 1) + 1;
  1946.                 denorm = 1;
  1947.                 }
  1948. #endif
  1949.         ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
  1950.         k = (int)ds;
  1951.         if (ds < 0. && ds != k)
  1952.                 k--;    /* want k = floor(ds) */
  1953.         k_check = 1;
  1954.         if (k >= 0 && k <= Ten_pmax) {
  1955.                 if (d < tens[k])
  1956.                         k--;
  1957.                 k_check = 0;
  1958.                 }
  1959.         j = bbits - i - 1;
  1960.         if (j >= 0) {
  1961.                 b2 = 0;
  1962.                 s2 = j;
  1963.                 }
  1964.         else {
  1965.                 b2 = -j;
  1966.                 s2 = 0;
  1967.                 }
  1968.         if (k >= 0) {
  1969.                 b5 = 0;
  1970.                 s5 = k;
  1971.                 s2 += k;
  1972.                 }
  1973.         else {
  1974.                 b2 -= k;
  1975.                 b5 = -k;
  1976.                 s5 = 0;
  1977.                 }
  1978.         if (mode < 0 || mode > 9)
  1979.                 mode = 0;
  1980.         try_quick = 1;
  1981.         if (mode > 5) {
  1982.                 mode -= 4;
  1983.                 try_quick = 0;
  1984.                 }
  1985.         leftright = 1;
  1986.         switch(mode) {
  1987.                 case 0:
  1988.                 case 1:
  1989.                         ilim = ilim1 = -1;
  1990.                         i = 18;
  1991.                         ndigits = 0;
  1992.                         break;
  1993.                 case 2:
  1994.                         leftright = 0;
  1995.                         /* no break */
  1996.                 case 4:
  1997.                         if (ndigits <= 0)
  1998.                                 ndigits = 1;
  1999.                         ilim = ilim1 = i = ndigits;
  2000.                         break;
  2001.                 case 3:
  2002.                         leftright = 0;
  2003.                         /* no break */
  2004.                 case 5:
  2005.                         i = ndigits + k + 1;
  2006.                         ilim = i;
  2007.                         ilim1 = i - 1;
  2008.                         if (i <= 0)
  2009.                                 i = 1;
  2010.                 }
  2011.         j = sizeof(unsigned long);
  2012.         for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i;
  2013.                 j <<= 1) result_k++;
  2014.         result = Balloc(result_k);
  2015.         s = s0 = (char *)result;
  2016.  
  2017.         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
  2018.  
  2019.                 /* Try to get by with floating-point arithmetic. */
  2020.  
  2021.                 i = 0;
  2022.                 d2 = d;
  2023.                 k0 = k;
  2024.                 ilim0 = ilim;
  2025.                 ieps = 2; /* conservative */
  2026.                 if (k > 0) {
  2027.                         ds = tens[k&0xf];
  2028.                         j = k >> 4;
  2029.                         if (j & Bletch) {
  2030.                                 /* prevent overflows */
  2031.                                 j &= Bletch - 1;
  2032.                                 d /= bigtens[n_bigtens-1];
  2033.                                 ieps++;
  2034.                                 }
  2035.                         for(; j; j >>= 1, i++)
  2036.                                 if (j & 1) {
  2037.                                         ieps++;
  2038.                                         ds *= bigtens[i];
  2039.                                         }
  2040.                         d /= ds;
  2041.                         }
  2042.                 else if (j1 = -k) {
  2043.                         d *= tens[j1 & 0xf];
  2044.                         for(j = j1 >> 4; j; j >>= 1, i++)
  2045.                                 if (j & 1) {
  2046.                                         ieps++;
  2047.                                         d *= bigtens[i];
  2048.                                         }
  2049.                         }
  2050.                 if (k_check && d < 1. && ilim > 0) {
  2051.                         if (ilim1 <= 0)
  2052.                                 goto fast_failed;
  2053.                         ilim = ilim1;
  2054.                         k--;
  2055.                         d *= 10.;
  2056.                         ieps++;
  2057.                         }
  2058.                 eps = ieps*d + 7.;
  2059.                 word0(eps) -= (P-1)*Exp_msk1;
  2060.                 if (ilim == 0) {
  2061.                         S = mhi = 0;
  2062.                         d -= 5.;
  2063.                         if (d > eps)
  2064.                                 goto one_digit;
  2065.                         if (d < -eps)
  2066.                                 goto no_digits;
  2067.                         goto fast_failed;
  2068.                         }
  2069. #ifndef No_leftright
  2070.                 if (leftright) {
  2071.                         /* Use Steele & White method of only
  2072.                          * generating digits needed.
  2073.                          */
  2074.                         eps = 0.5/tens[ilim-1] - eps;
  2075.                         for(i = 0;;) {
  2076.                                 L = (long)d;
  2077.                                 d -= L;
  2078.                                 *s++ = '0' + (int)L;
  2079.                                 if (d < eps)
  2080.                                         goto ret1;
  2081.                                 if (1. - d < eps)
  2082.                                         goto bump_up;
  2083.                                 if (++i >= ilim)
  2084.                                         break;
  2085.                                 eps *= 10.;
  2086.                                 d *= 10.;
  2087.                                 }
  2088.                         }
  2089.                 else {
  2090. #endif
  2091.                         /* Generate ilim digits, then fix them up. */
  2092.                         eps *= tens[ilim-1];
  2093.                         for(i = 1;; i++, d *= 10.) {
  2094.                                 L = (long)d;
  2095.                                 d -= L;
  2096.                                 *s++ = '0' + (int)L;
  2097.                                 if (i == ilim) {
  2098.                                         if (d > 0.5 + eps)
  2099.                                                 goto bump_up;
  2100.                                         else if (d < 0.5 - eps) {
  2101.                                                 while(*--s == '0');
  2102.                                                 s++;
  2103.                                                 goto ret1;
  2104.                                                 }
  2105.                                         break;
  2106.                                         }
  2107.                                 }
  2108. #ifndef No_leftright
  2109.                         }
  2110. #endif
  2111.  fast_failed:
  2112.                 s = s0;
  2113.                 d = d2;
  2114.                 k = k0;
  2115.                 ilim = ilim0;
  2116.                 }
  2117.  
  2118.         /* Do we have a "small" integer? */
  2119.  
  2120.         if (be >= 0 && k <= Int_max) {
  2121.                 /* Yes. */
  2122.                 ds = tens[k];
  2123.                 if (ndigits < 0 && ilim <= 0) {
  2124.                         S = mhi = 0;
  2125.                         if (ilim < 0 || d <= 5*ds)
  2126.                                 goto no_digits;
  2127.                         goto one_digit;
  2128.                         }
  2129.                 for(i = 1;; i++) {
  2130.                         L = (long)(d / ds);
  2131.                         d -= L*ds;
  2132. #ifdef Check_FLT_ROUNDS
  2133.                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
  2134.                         if (d < 0) {
  2135.                                 L--;
  2136.                                 d += ds;
  2137.                                 }
  2138. #endif
  2139.                         *s++ = '0' + (int)L;
  2140.                         if (i == ilim) {
  2141.                                 d += d;
  2142.                                 if (d > ds || d == ds && L & 1) {
  2143.  bump_up:
  2144.                                         while(*--s == '9')
  2145.                                                 if (s == s0) {
  2146.                                                         k++;
  2147.                                                         *s = '0';
  2148.                                                         break;
  2149.                                                         }
  2150.                                         ++*s++;
  2151.                                         }
  2152.                                 break;
  2153.                                 }
  2154.                         if (!(d *= 10.))
  2155.                                 break;
  2156.                         }
  2157.                 goto ret1;
  2158.                 }
  2159.  
  2160.         m2 = b2;
  2161.         m5 = b5;
  2162.         mhi = mlo = 0;
  2163.         if (leftright) {
  2164.                 if (mode < 2) {
  2165.                         i =
  2166. #ifndef Sudden_Underflow
  2167.                                 denorm ? be + (Bias + (P-1) - 1 + 1) :
  2168. #endif
  2169. #ifdef IBM
  2170.                                 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
  2171. #else
  2172.                                 1 + P - bbits;
  2173. #endif
  2174.                         }
  2175.                 else {
  2176.                         j = ilim - 1;
  2177.                         if (m5 >= j)
  2178.                                 m5 -= j;
  2179.                         else {
  2180.                                 s5 += j -= m5;
  2181.                                 b5 += j;
  2182.                                 m5 = 0;
  2183.                                 }
  2184.                         if ((i = ilim) < 0) {
  2185.                                 m2 -= i;
  2186.                                 i = 0;
  2187.                                 }
  2188.                         }
  2189.                 b2 += i;
  2190.                 s2 += i;
  2191.                 mhi = i2b(1);
  2192.                 }
  2193.         if (m2 > 0 && s2 > 0) {
  2194.                 i = m2 < s2 ? m2 : s2;
  2195.                 b2 -= i;
  2196.                 m2 -= i;
  2197.                 s2 -= i;
  2198.                 }
  2199.         if (b5 > 0) {
  2200.                 if (leftright) {
  2201.                         if (m5 > 0) {
  2202.                                 mhi = pow5mult(mhi, m5);
  2203.                                 b1 = mult(mhi, b);
  2204.                                 Bfree(b);
  2205.                                 b = b1;
  2206.                                 }
  2207.                         if (j = b5 - m5)
  2208.                                 b = pow5mult(b, j);
  2209.                         }
  2210.                 else
  2211.                         b = pow5mult(b, b5);
  2212.                 }
  2213.         S = i2b(1);
  2214.         if (s5 > 0)
  2215.                 S = pow5mult(S, s5);
  2216.  
  2217.         /* Check for special case that d is a normalized power of 2. */
  2218.  
  2219.         if (mode < 2) {
  2220.                 if (!word1(d) && !(word0(d) & Bndry_mask)
  2221. #ifndef Sudden_Underflow
  2222.                  && word0(d) & Exp_mask
  2223. #endif
  2224.                                 ) {
  2225.                         /* The special case */
  2226.                         b2 += Log2P;
  2227.                         s2 += Log2P;
  2228.                         spec_case = 1;
  2229.                         }
  2230.                 else
  2231.                         spec_case = 0;
  2232.                 }
  2233.  
  2234.         /* Arrange for convenient computation of quotients:
  2235.          * shift left if necessary so divisor has 4 leading 0 bits.
  2236.          *
  2237.          * Perhaps we should just compute leading 28 bits of S once
  2238.          * and for all and pass them and a shift to quorem, so it
  2239.          * can do shifts and ors to compute the numerator for q.
  2240.          */
  2241. #ifdef Pack_32
  2242.         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
  2243.                 i = 32 - i;
  2244. #else
  2245.         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
  2246.                 i = 16 - i;
  2247. #endif
  2248.         if (i > 4) {
  2249.                 i -= 4;
  2250.                 b2 += i;
  2251.                 m2 += i;
  2252.                 s2 += i;
  2253.                 }
  2254.         else if (i < 4) {
  2255.                 i += 28;
  2256.                 b2 += i;
  2257.                 m2 += i;
  2258.                 s2 += i;
  2259.                 }
  2260.         if (b2 > 0)
  2261.                 b = lshift(b, b2);
  2262.         if (s2 > 0)
  2263.                 S = lshift(S, s2);
  2264.         if (k_check) {
  2265.                 if (cmp(b,S) < 0) {
  2266.                         k--;
  2267.                         b = multadd(b, 10, 0);  /* we botched the k estimate */
  2268.                         if (leftright)
  2269.                                 mhi = multadd(mhi, 10, 0);
  2270.                         ilim = ilim1;
  2271.                         }
  2272.                 }
  2273.         if (ilim <= 0 && mode > 2) {
  2274.                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
  2275.                         /* no digits, fcvt style */
  2276.  no_digits:
  2277.                         k = -1 - ndigits;
  2278.                         goto ret;
  2279.                         }
  2280.  one_digit:
  2281.                 *s++ = '1';
  2282.                 k++;
  2283.                 goto ret;
  2284.                 }
  2285.         if (leftright) {
  2286.                 if (m2 > 0)
  2287.                         mhi = lshift(mhi, m2);
  2288.  
  2289.                 /* Compute mlo -- check for special case
  2290.                  * that d is a normalized power of 2.
  2291.                  */
  2292.  
  2293.                 mlo = mhi;
  2294.                 if (spec_case) {
  2295.                         mhi = Balloc(mhi->k);
  2296.                         Bcopy(mhi, mlo);
  2297.                         mhi = lshift(mhi, Log2P);
  2298.                         }
  2299.  
  2300.                 for(i = 1;;i++) {
  2301.                         dig = quorem(b,S) + '0';
  2302.                         /* Do we yet have the shortest decimal string
  2303.                          * that will round to d?
  2304.                          */
  2305.                         j = cmp(b, mlo);
  2306.                         delta = diff(S, mhi);
  2307.                         j1 = delta->sign ? 1 : cmp(b, delta);
  2308.                         Bfree(delta);
  2309. #ifndef ROUND_BIASED
  2310.                         if (j1 == 0 && !mode && !(word1(d) & 1)) {
  2311.                                 if (dig == '9')
  2312.                                         goto round_9_up;
  2313.                                 if (j > 0)
  2314.                                         dig++;
  2315.                                 *s++ = dig;
  2316.                                 goto ret;
  2317.                                 }
  2318. #endif
  2319.                         if (j < 0 || j == 0 && !mode
  2320. #ifndef ROUND_BIASED
  2321.                                                         && !(word1(d) & 1)
  2322. #endif
  2323.                                         ) {
  2324.                                 if (j1 > 0) {
  2325.                                         b = lshift(b, 1);
  2326.                                         j1 = cmp(b, S);
  2327.                                         if ((j1 > 0 || j1 == 0 && dig & 1)
  2328.                                         && dig++ == '9')
  2329.                                                 goto round_9_up;
  2330.                                         }
  2331.                                 *s++ = dig;
  2332.                                 goto ret;
  2333.                                 }
  2334.                         if (j1 > 0) {
  2335.                                 if (dig == '9') { /* possible if i == 1 */
  2336.  round_9_up:
  2337.                                         *s++ = '9';
  2338.                                         goto roundoff;
  2339.                                         }
  2340.                                 *s++ = dig + 1;
  2341.                                 goto ret;
  2342.                                 }
  2343.                         *s++ = dig;
  2344.                         if (i == ilim)
  2345.                                 break;
  2346.                         b = multadd(b, 10, 0);
  2347.                         if (mlo == mhi)
  2348.                                 mlo = mhi = multadd(mhi, 10, 0);
  2349.                         else {
  2350.                                 mlo = multadd(mlo, 10, 0);
  2351.                                 mhi = multadd(mhi, 10, 0);
  2352.                                 }
  2353.                         }
  2354.                 }
  2355.         else
  2356.                 for(i = 1;; i++) {
  2357.                         *s++ = dig = quorem(b,S) + '0';
  2358.                         if (i >= ilim)
  2359.                                 break;
  2360.                         b = multadd(b, 10, 0);
  2361.                         }
  2362.  
  2363.         /* Round off last digit */
  2364.  
  2365.         b = lshift(b, 1);
  2366.         j = cmp(b, S);
  2367.         if (j > 0 || j == 0 && dig & 1) {
  2368.  roundoff:
  2369.                 while(*--s == '9')
  2370.                         if (s == s0) {
  2371.                                 k++;
  2372.                                 *s++ = '1';
  2373.                                 goto ret;
  2374.                                 }
  2375.                 ++*s++;
  2376.                 }
  2377.         else {
  2378.                 while(*--s == '0');
  2379.                 s++;
  2380.                 }
  2381.  ret:
  2382.         Bfree(S);
  2383.         if (mhi) {
  2384.                 if (mlo && mlo != mhi)
  2385.                         Bfree(mlo);
  2386.                 Bfree(mhi);
  2387.                 }
  2388.  ret1:
  2389.         Bfree(b);
  2390.         *s = 0;
  2391.         *decpt = k + 1;
  2392.         if (rve)
  2393.                 *rve = s;
  2394.         return s0;
  2395.         }
  2396.  
  2397. #ifdef atarist
  2398. double atof(const char *s)
  2399. {
  2400.     return strtod(s, (char **)NULL);
  2401. }
  2402. #endif
  2403.  
  2404. #endif /* USE_DTOA */
  2405.